
#include <math.h>

#include "globals.h"
#include "geometry.h"
#include "distance.h"
#include "sgidrawi.h"

void
perpendicularLine (Pt *point, Line *line, Line *perp)
{
    perp->w = (point->x)*(line->y) - (point->y)*(line->x);
    perp->x = (point->w)*(line->y)*-1;
    perp->y = (point->w)*(line->x);
}

void
perpRightTurnLine (Pt *point, Line *line, Line *perp)
{
	/*	Give the point a positive w value. */

	if (point->w < 0) {
		point->w *= -1;
		point->x *= -1;
		point->y *= -1;
	}

	/*	Get the perpendicular line, right turn from "line". */

    perp->w = (point->y)*(line->x) - (point->x)*(line->y);
    perp->x = (point->w)*(line->y);
    perp->y = (point->w)*(line->x)*-1;

	/*	Lastly, translate the line so that it passes through the point.  */

	/*	The above set-up also moved the line to the point. */

}

void
ptNearestLine(Pt *point, Line *line, Pt *nearest)
{
	Line	perp;

	perpendicularLine (point, line, &perp);
    IntersectionOf(line, &perp, nearest);
	RoundPoint(*nearest);
}

boolean
pastSegmentBounds(Pt *point, Pt *start, Pt *end, int *closerEnd)
{
	Real		dist0, dist1;
	boolean		retVal;
	int			indicator1, indicator2;

    squarePointDistance(point, start, &dist0);
    squarePointDistance(point, end, &dist1);

    indicator1 = SignOfaLessb(XValPtr(point), XValPtr(start)) *
          SignOfaLessb(XValPtr(point), XValPtr(end));
    indicator2 = SignOfaLessb(YValPtr(point), YValPtr(start)) *
          SignOfaLessb(YValPtr(point), YValPtr(end));

    if ( (indicator1 > 0) || ((indicator1 == 0) && (indicator2 > 0)) ) {

		retVal = True;

        /*  Both direction vectors go in the same direction so
            the "nearest" point does not lie between the start and end. */

        if (Less(dist0, dist1)) 
			*closerEnd = 0;
        else 
			*closerEnd = 1;
    }
	else
		retVal = False;

	return (retVal);
}

void
ptNearestSegment(Pt *point, Pt *start, Pt *end, Pt *nearest)
{
	Line	fullLine;
	int		nearEnd;

	LineOf(start,end,&fullLine);
	ptNearestLine(point, &fullLine, nearest);

	/*	Now, check to see if we have gone beyond the endpoints. */

	if (pastSegmentBounds(nearest, start, end, &nearEnd) == True) {

        /*  Both direction vectors go in the same direction so
            the "nearest" point does not lie between the start and end. */

        if (nearEnd == 0) {
			PointAssignPtr(nearest, start);
		}
        else {
			PointAssignPtr(nearest, end);
		}
    }
}

void
createParabola(Pt *focus, Line *directrix, Parabola *para)
{
	Pt		nearPt, transfocus;
	double	twoFocalLength;
	Line	perp;

#ifdef DEBUG_GEOM
	printf ("geom: find parabola of focus %f %f %f and line %f %f %f\n",
		focus->w, focus->x, focus->y,
		directrix->w, directrix->x, directrix->y);
#endif

	if (OnTheLine(directrix, focus) == False) {
		para->scale = 1.0;
		para->offset = 0.0;
		para->sqCoeff = 1;
		ptNearestLine (focus, directrix, &nearPt);
		pointDistance(focus, &nearPt, &twoFocalLength);

		MidPointPtr(focus, &nearPt, &(para->centre));
		RoundPoint(para->centre);
		PointAssignPtr(&(para->focus), focus);
		PointAssignPtr(&(para->directrix), directrix);
		para->directrix.w = 0.0;
		RoundPoint(para->focus);

		para->focalLength = (Real)(twoFocalLength / 2);
			/*	For the sine of the angle, we shuld multiply by -1 to get
				the cosine (due to the 90 degree angle difference), but
				I want the negative of the sine, so two multiplications by
				-1 cancel each other out. */
		para->dirSin = 
			(Real) MyAbs(((XValPtr(focus) - XVal(nearPt)) /twoFocalLength));
		para->dirCos = 
			(Real) MyAbs(((YValPtr(focus) - YVal(nearPt)) /twoFocalLength));

		transfocus.w = focus->w * nearPt.w;
		transfocus.x = focus->x * nearPt.w - focus->w * nearPt.x;
		transfocus.y = focus->y * nearPt.w - focus->w * nearPt.y;

		if ((XVal(transfocus) >= 0) && (YVal(transfocus) > 0)) {
		}
		else if ((XVal(transfocus) < 0) && (YVal(transfocus) > 0)) {
			para->dirSin *= -1;
		}
		else if ((XVal(transfocus) < 0) && (YVal(transfocus) <= 0)) {
			para->dirSin *= -1;
			para->dirCos *= -1;
		}
		else {
			para->dirCos *= -1;
		}
	}
	else {
		para->scale = 1.0;
		para->offset = 0.0;
		para->sqCoeff = 0;

		PointAssignPtr(&(para->focus), focus);
		PointAssignPtr(&(para->directrix), directrix);
		para->directrix.w = 0.0;
		PointAssignPtr(&(para->centre), focus);
		RoundPoint(para->centre);
		RoundPoint(para->focus);

		/*	Encode the perpendicular line in the focal length and
			direction parameters. */

		perpendicularLine(focus, directrix, &perp);

		para->focalLength = (Real) perp.w;
		para->dirSin = (Real) perp.x;
		para->dirCos = (Real) perp.y;
	}
	PointAssign(para->scalePt[0], para->centre);
	PointAssign(para->scalePt[1], para->centre);
}

Real
parabolaParam(Parabola *para, Pt *pt)
{
	return((Real) ((XValPtr(pt) - XVal(para->centre)) * para->dirCos -
			(YValPtr(pt) - YVal(para->centre)) * para->dirSin));
}

void
parabolaPoint(Parabola *para, Real param, Pt *pt)
{
	pt->w = 1.0;
	pt->x = param*(para->dirCos + param * para->dirSin / (4*para->focalLength))
			+ XVal(para->centre);
	pt->y = param*(param * para->dirCos / (4*para->focalLength) - para->dirSin)
			+ YVal(para->centre);
}

void
intersectLineParabola(Line *line, Parabola *para, Pt *point1, Pt *point2)
{
	Real	a,b,c, soln1, soln2, discriminant, rootDisc;
	Real	squared; 
	Line	perp, trans;

	if (para->sqCoeff == 0) {
		/*	We have a degenerate parabola.  The intersection is easy. */
		perp.w = para->focalLength;
		perp.x = para->dirSin;
		perp.y = para->dirCos;

		IntersectionOf(line, &perp, point1);
		PointAssignPtr(point2, point1);	
	}
	else {
		/*	The intersection is characterised by a quadratic equation
			constraint on the parameter for the parabola.  Start by 
			setting out the coefficients for this quadratic equation:
				a*t^2 + b*t + c = 0			*/

		a = (line->x * para->dirSin + line->y * para->dirCos) / 
				(4 * para->focalLength);
		b = line->x * para->dirCos - line->y * para->dirSin;
		c = line->w + (line->x * XVal(para->centre)) + 
				(line->y * YVal(para->centre));
		PointAssignPtr(&trans, line);
		trans.w = 0;

		if ((!Zero(a)) && (!Zero(DotProd(&trans, &(para->directrix))))) {
			/*	Find the discriminant of the equation.  We assume that some
				intersection will occur, so we use a zero discriminant if
				it turns out to be negative. */

			discriminant = b*b -4*a*c;
			discriminant = (Zero(discriminant) ? Max(discriminant, 0) :
								discriminant);

			/*	mjm's hack for now. */
			discriminant = Max(discriminant, 0);

			if (discriminant >= 0) {
				rootDisc = sqrt((double)discriminant);
				soln1 = (b * -1 + rootDisc) / (a * 2);
				soln2 = (b * -1 - rootDisc) / (a * 2);

#ifdef DEBUG_GEOM
				printf ("geo: coefficients are a %2.2f, b %2.2f, c %2.2f\n", 
					a, b, c);
				printf ("geo: discriminant is %2.2f with solutions %f and %f\n",
					discriminant, soln1, soln2);
#endif

				/*	Find the solution points. */

				parabolaPoint(para, soln1, point1);
				parabolaPoint(para, soln2, point2);

#ifdef DEBUG_GEOM
				printf ("geo: candidate points are %2.2f %2.2f "
					"and %2.2f %2.2f\n", XValPtr(point1), YValPtr(point1), 
					XValPtr(point2), YValPtr(point2));
#endif
			}
			else {
				/* mjm to fill in later. */
#if 0
				printf ("Oh no!  Parabola and line do not intersect!\n");
				printf ("Line is %f %f %f\n", line->w, line->x, line->y);
				printf ("Parabola is focus %f %f, centre %f %f\n",
					XVal(para->focus), YVal(para->focus),
					XVal(para->centre), YVal(para->centre));
				printf ("  focal %f, coeff %d, scale %f, offset %f\n",
					para->focalLength, para->sqCoeff, para->scale,
					para->offset);
				printf ("  Sin %f, cos %f\n", para->dirSin, para->dirCos);
				printf ("  Discriminant was %f\n", discriminant);
#endif
			}
		}
		else {
			soln1 = c * -1 / b;

			squared = soln1 * soln1;
			point1->w = 1.0;
			point1->x = soln1* para->dirCos + squared * para->dirSin /
						(para->focalLength * 4) + XVal(para->centre);
			point1->y = squared * para->dirCos / (para->focalLength * 4) - 
						soln1 * para->dirSin + YVal(para->centre);

			PointAssignPtr(point2, point1);

#ifdef DEBUG_GEOM
			printf ("geo: single parameter solution at %f, point %f %f %f\n",
				soln1, point1->w, point1->x, point1->y);
#endif

		}
	}
}

/*	Adjust the points on a parabola so that they are just within the 
	display window.  This way, we don't waste time drawing stuff we
	don't see.  For now, the routine takes advantage of the fact that
	it should only be called with the sweepline so I know the 
	orientation of the parabolas. */

void
screenTrim(Parabola *para, Pt *pt1, Pt *pt2)
{
	Line	top, bottom, left, right;
	Real	topParm, bottomParm, leftParm, rightParm;
	Pt		topPt, bottomPt, leftPt, rightPt, dummy;

#if 0
printf ("Reorient in:  %f %f and %f %f\n",  XValPtr(pt1), YValPtr(pt1),
XValPtr(pt2), YValPtr(pt2));
	top.w = -125.0; top.x = 0.0; top.y = 1.0;
	bottom.w = 25.0; bottom.x = 0.0; bottom.y = 1.0;
	right.w = -125.0; right.x = 1.0; right.y = 0.0;
	left.w = 25.0; left.x = 1.0; left.y = 0.0;

	intersectLineParabola(&top, para, &topPt, &dummy);
	intersectLineParabola(&bottom, para, &bottomPt, &dummy);
	intersectLineParabola(&left, para, &leftPt, &dummy);
	intersectLineParabola(&right, para, &rightPt, &dummy);

	/*	First, bring it into the acceptable Y range. */

	if (YValPtr(pt1) > YVal(topPt)) {
		PointAssign(*pt1, topPt);
	}

	if (YValPtr(pt2) > YVal(topPt)) {
		PointAssign(*pt2, topPt);
	}

	if (YValPtr(pt1) < YVal(bottomPt)) {
		PointAssign(*pt1, bottomPt);
	}

	if (YValPtr(pt2) < YVal(bottomPt)) {
		PointAssign(*pt2, bottomPt);
	}

	/*	Then, bring it into the acceptable X range. */

	if (XValPtr(pt1) > XVal(rightPt)) {
		PointAssign(*pt1, rightPt);
	}

	if (XValPtr(pt2) > XVal(rightPt)) {
		PointAssign(*pt2, rightPt);
	}

	if (XValPtr(pt1) < XVal(leftPt)) {
		PointAssign(*pt1, leftPt);
	}

	if (XValPtr(pt2) < XVal(leftPt)) {
		PointAssign(*pt2, leftPt);
	}
printf ("Reorient out:  %f %f and %f %f\n",  XValPtr(pt1), YValPtr(pt1),
XValPtr(pt2), YValPtr(pt2));
#endif
}

void
reorientParabola (Parabola *para, Pt *point1, Pt *point2, boolean trimToScreen)
{
	Real	x1, x2;
	Real	small, big;
	Pt		newPt1, newPt2;

	if (para->sqCoeff != 0) {
		/*	First see if the points should be tailored to the screen itself.  */

		PointAssign(newPt1, *point1);
		PointAssign(newPt2, *point2);
		if (trimToScreen == True) {
			screenTrim(para, &newPt1, &newPt2);
		}

		/*	Reorient the two points to standard form to find the first
			one along the parameterization of the parabola. */

		x1 = parabolaParam(para, point1);
		x2 = parabolaParam(para, point2);

		if (Less(x1, x2)) {
			small = x1;
			big = x2;
			PointAssign(para->scalePt[0], newPt1);
			PointAssign(para->scalePt[1], newPt2);
		}
		else {
			small = x2;
			big = x1;
			PointAssign(para->scalePt[0], newPt2);
			PointAssign(para->scalePt[1], newPt1);
		}

		/*	We can now set the scale factors. */

		para->scale = big - small;
		para->offset = small;
	}
	else {
		PointAssign(para->scalePt[0], *point1);
		PointAssign(para->scalePt[1], *point2);
	}
}

void
drawParabola(int colour, Parabola *para)
{
#ifdef Graphics
	int		i;
	float	param;
	int		divisions;
	Pt		pt0, pt1;

	if (para->sqCoeff != 0) {
		if (para->focalLength > 0.0001) {
			param = para->offset;
			PointAssign(pt1, para->scalePt[0]);
			divisions = (para->scale <  50 ? 10 : 
						(para->scale < 300 ? 15 :
						(para->scale < 600 ? 20 : 
											 25)));
			for (i = 1; i <= divisions; i++) {
				PointAssign(pt0, pt1);
				if (i < divisions) {
					param = para->scale * i / divisions + para->offset;
					parabolaPoint(para, param, &pt1);
				}
				else {
					param = para->scale + para->offset;
					PointAssign(pt1, para->scalePt[1]);
				}
				drawline(colour, XVal(pt0), YVal(pt0), XVal(pt1), YVal(pt1));
			}
		}
		else {
			if (para->offset * (para->offset + para->scale) < 0) {
				drawline(colour,XVal(para->scalePt[0]), YVal(para->scalePt[0]), 
					XVal(para->centre), YVal(para->centre));
				drawline(colour, XVal(para->centre), YVal(para->centre),
					XVal(para->scalePt[1]), YVal(para->scalePt[1]));
			}
			else {
				drawline(colour,XVal(para->scalePt[0]), YVal(para->scalePt[0]), 
					XVal(para->scalePt[1]), YVal(para->scalePt[1]));
			}
		}
	}
	else {
		if (SamePoint(&(para->scalePt[0]), &(para->scalePt[1]))) {
			drawline(colour, XVal(para->scalePt[0]), YVal(para->scalePt[0]), 
					XVal(para->centre), YVal(para->centre));
		}
		else {
			drawline(colour, XVal(para->scalePt[0]), YVal(para->scalePt[0]), 
					XVal(para->scalePt[1]), YVal(para->scalePt[1]));
		}
	}
#endif
}

